package fem2.element;

import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import junit.framework.Assert;

import org.junit.Test;

import fem2.Constraint;
import fem2.Debugger;
import fem2.Element;
import fem2.Face;
import fem2.Mesh;
import fem2.Model;
import fem2.Node;
import fem2.analysis.Analysis;
import fem2.analysis.LinearStaticAnalysis;
import fem2.enu.EchoLevelType;
import fem2.enu.SolverType;

public class PoissonElementTest {

	public Element createElement() {
		Node n1 = new Node(0.0, 0.0);
		Node n2 = new Node(1.0, 0.0);
		Node n3 = new Node(1.0, 1.0);
		Node n4 = new Node(0.0, 1.0);
		Face f = new Face(n1, n2, n3, n4);

		f.setNumberOfFields(1);
		n1.setNumberOfFields(1);
		n2.setNumberOfFields(1);
		n3.setNumberOfFields(1);
		n4.setNumberOfFields(1);

		return new PoissonElement(f, 1.0);
	}

	@Test
	public void test1() {
		Element e = createElement();

		Face f = (Face) e.getMeshPart();
		int nDofs = 0;
		for (int i = 0; i < 4; i++) {
			nDofs = f.getNodes()[i].enumerateDofs(nDofs);
		}
		f.enumerateDofs();

		e.Initialize();

		IMatrix kte = new Array2DMatrix(4, 4);
		e.computeKte(kte);

		Debugger.watch(kte);

		/*
		 * column 1
		 */
		double[] uh1 = new double[] { 1e-3, 0, 0, 0 };
		for (int i = 0; i < 4; i++) {
			f.getNodes()[i].setUHatGlobal(0.0, uh1);
		}
		f.setUHatGlobal(0.0, uh1);

		e.FinalizeNonlinearIteration();

		IVector rie = new ArrayVector(4);
		e.computeRie(rie);
		BLAM.scale(1 / 1e-3, rie);

		Debugger.watch(rie);

		/*
		 * column 2
		 */
		double[] uh2 = new double[] { 0, 1e-3, 0, 0 };
		for (int i = 0; i < 4; i++) {
			f.getNodes()[i].setUHatGlobal(0.0, uh2);
		}
		f.setUHatGlobal(0.0, uh2);

		e.FinalizeNonlinearIteration();

		BLAM.scale(0.0, rie);
		e.computeRie(rie);
		BLAM.scale(1 / 1e-3, rie);

		Debugger.watch(rie);

		/*
		 * column 3
		 */
		double[] uh3 = new double[] { 0, 0, 1e-3, 0 };
		for (int i = 0; i < 4; i++) {
			f.getNodes()[i].setUHatGlobal(0.0, uh3);
		}
		f.setUHatGlobal(0.0, uh3);

		e.FinalizeNonlinearIteration();

		BLAM.scale(0.0, rie);
		e.computeRie(rie);
		BLAM.scale(1 / 1e-3, rie);

		Debugger.watch(rie);

		/*
		 * column 4
		 */
		double[] uh4 = new double[] { 0, 0, 0, 1e-3 };
		for (int i = 0; i < 4; i++) {
			f.getNodes()[i].setUHatGlobal(0.0, uh4);
		}
		f.setUHatGlobal(0.0, uh4);

		e.FinalizeNonlinearIteration();

		BLAM.scale(0.0, rie);
		e.computeRie(rie);
		BLAM.scale(1 / 1e-3, rie);

		Debugger.watch(rie);
	}

	public Model createModel() {
		Mesh mesh = new Mesh();

		Node n1 = mesh.addNode(new Node(0, 0));
		Node n2 = mesh.addNode(new Node(0.5, 0));
		Node n3 = mesh.addNode(new Node(1, 0));

		Node n4 = mesh.addNode(new Node(0, 0.5));
		Node n5 = mesh.addNode(new Node(0.5, 0.5));
		Node n6 = mesh.addNode(new Node(1, 0.5));

		Node n7 = mesh.addNode(new Node(0, 1));
		Node n8 = mesh.addNode(new Node(0.5, 1));
		Node n9 = mesh.addNode(new Node(1, 1));

		Face f1 = mesh.addFace(new Face(n1, n2, n5, n4));
		Face f2 = mesh.addFace(new Face(n2, n3, n6, n5));
		Face f3 = mesh.addFace(new Face(n5, n6, n9, n8));
		Face f4 = mesh.addFace(new Face(n4, n5, n8, n7));

		Model m = new Model(mesh);

		m.addElement(new PoissonElement(f1, 0.3));
		m.addElement(new PoissonElement(f2, 0.3));
		m.addElement(new PoissonElement(f3, 0.3));
		m.addElement(new PoissonElement(f4, 0.3));

		Constraint c1 = new Constraint(false);
		c1.setValue(0, 0.0);

		n1.setConstraint(c1);
		n4.setConstraint(c1);
		n7.setConstraint(c1);

		Constraint c2 = new Constraint(false);
		c2.setValue(0, 5.0);

		n3.setConstraint(c2);
		n6.setConstraint(c2);
		n9.setConstraint(c2);

		return m;
	}

	@Test
	public void testModel1() {
		Model m = createModel();
		Analysis an = new LinearStaticAnalysis(m);
		an.getStrategy().setSolverType(SolverType.BAITSCH_SOLVER_GENERIC);
		an.getStrategy()
				.setEchoLevel(EchoLevelType.OUTPUT, EchoLevelType.AUX1, EchoLevelType.DEBUG);

		an.run();

		/*
		 * post process
		 */
		Node n2 = m.getMesh().getNode(1);
		Node n5 = m.getMesh().getNode(4);
		Node n8 = m.getMesh().getNode(7);

		System.out.println("T at node 2 = " + n2.getUHatLocal()[0]);
		System.out.println("T at node 5 = " + n5.getUHatLocal()[0]);
		System.out.println("T at node 8 = " + n8.getUHatLocal()[0]);

		Assert.assertEquals(2.5, n2.getUHatLocal()[0], 1e-10);
		Assert.assertEquals(2.5, n5.getUHatLocal()[0], 1e-10);
		Assert.assertEquals(2.5, n8.getUHatLocal()[0], 1e-10);
	}
}
